#ifndef _SYSTEM_M_HF_H_
#define _SYSTEM_M_HF_H_
#include<eigen3/Eigen/Eigen>//matrix library

#include"hfsolver_m.h"

using Eigen::Matrix;
using Eigen::Dynamic;
using Eigen::SelfAdjointEigenSolver;

 
/// system under HF basis
class System_TzM_HF:public System_TzM<HF_Orbital_TzM,HF_Orbital_TzM::DataType>
{
 public:
  typedef System_TzM_HO System_TzType;
  typedef HFSolver_M HFSolverType;
  typedef typename System_TzType::DataType DataType;

  typedef Matrix<DataType,Dynamic,Dynamic> Mat;
  typedef Matrix<DataType,Dynamic,1> Vec;
  
  typedef HF_Orbital_TzM OrbitalType;
  const static int MatEledim;//< dimension of TwoBodyMatEle.
  typedef TwoBodyMatElement<DataType> TwoBodyMatEle;
  typedef Matrix<TwoBodyMatEle,Dynamic,Dynamic> Mat2B;

  System_TzM_HF(){}
 System_TzM_HF(System_TzType * Ptr):pSystem_Tz(Ptr),hfsolver(Ptr),hbar_omega(Ptr->hbar_omega){}
 System_TzM_HF(System_TzType * Ptr,int _A,int _Z):pSystem_Tz(Ptr),hfsolver(Ptr,_A,_Z),A(_A),Z(_Z),hbar_omega(Ptr->hbar_omega){}

  void setAZ(int _A,int _Z)
  {
    A=_A;
    Z=_Z;
    hfsolver.setAZ(_A,_Z);
  }

  void setupOrbitals();
  void setupTwoBodyMat();

  void setup(System_TzType * Ptr,int _A,int _Z)
  {
    pSystem_Tz=Ptr;
    hbar_omega=Ptr->hbar_omega;
    A=_A;
    Z=_Z;
    hfsolver.setup(Ptr,_A,_Z);
    setup();
  }
  void setup()
  {
    setupOrbitals();
    setupTwoBodyStates();
    //setupTwoBodyMat();
    beta=0.;
  }
  ///interface for method beyond HF.
  DataType get1B(int bra,int ket) const;  
  DataType get2B(int a,int b,int c,int d) const
  {
    TwoBodyMatEle MatEle;
    get2BmatOnNeed(a,b,c,d,MatEle);
    //set2BmatAt(a,b,c,d,MatEle);
    if(MatEle.empty()) return 0;
    return MatEle[0]+(beta-1)*MatEle[3]*hbar_omega/A-beta*MatEle[2]*hbar_omega/A;    
  }
  DataType get0B() const
  {
    // return 0;
    return -hbar_omega*1.5*beta;    
  }

  ///interface for nocore shell model
  DataType get1B_p2(int bra,int ket) const;
  double get1B_p2_plus_r2(int bra,int ket) const;
  
  void get2BmatOnNeed(int a,int b,int c,int d,TwoBodyMatEle & MatEle) const;

  void printOrbitals(int num=-1) const
  {
    if(num==-1)
      num=Orbitals.size();
    cout<<"num"<<"\t"<<"par"<<"\t"<<"mm"<<"\t"<<"tz"<<"\t"<<"e"<<"\t"<<"OccuP"<<endl;
    for(int i=0;i<num;i++)
      {
      	cout<<i<<"\t"<<Orbitals[i].par<<"\t"<<Orbitals[i].mm<<"\t"<<Orbitals[i].tz<<"\t"<<Orbitals[i].e<<"\t"<<1-Orbitals[i].state<<endl;
      }
  }

  
  ///maybe usefull operator
  ///only radial part
  DataType rL(int bra,int ket,int L) const;
  /* DataType rms() const */
  /* { */
  /*   DataType val=0; */
  /*   for(int i=0;i<=FermiSurface;i++) */
  /*     { */
  /* 	val+=rL(i,i,2); */
  /*     } */
  /*   return val/A; */
  /* } */


  DataType jplus(int a,int b) const;
  DataType jminus(int a,int b) const;
  
  int A,Z;//for HF. basis generating.


  int maxProtons;
  int maxNeutrons;
  //  int FermiSurface;
  DataType HFenergy;
  System_TzType * pSystem_Tz;
  HFSolverType hfsolver;

  double hbar_omega;
  double beta;
 private:
};



#endif
